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ABSTRACT 



Aims. The goal of this paper is twofold. In the first place we complete the set of diagnostic tools for synchrotron emitting sources 
presented by Del Zanna et al. (Astron. Astrophys. 453, 621, 2006] with the computation of inverse Compton radiation from the same 
relativistic particles. Moreover we investigate, for the first time, the gamma-ray emission properties of Pulsar Wind Nebulae in the 

(-H ' light of the axisymmetricjef-toriii scenario. 

Oj. Methods. The proposed method consists in evolving the relativistic MHD equations and the maximum energy of the emitting particles 

I . including adiabatic and synchrotron losses along streamlines. The particle energy distribution function is split in two components: 

O ' one corresponding to the radio emitting electrons interpreted as a relic population born at the outburst of the supernova and the other 

associated with the wind population continuously accelerated at the termination shock and emitting up to the gamma-ray band. The 

C/5 , inverse Compton emissivity is calculated using the general Klein-Nishina differential cross-section and three different photon targets 

Cu . for the relativistic particles are considered: the nebular synchrotron photons, photons associated with the far-infrared thermal excess 

and the cosmic microwave background. 

Results. When the method is applied to the simulations that better reproduce the optical and X-ray morphology of the Crab Nebula, 
the overall synchrotron spectrum can only be fitted assuming an excess of injected particles and a steeper power law {E"-'') with 
respect to previous models. The resulting TeV emission has then the correct shape but is in excess of the data. This is related to the 
magnetic field structure in the nebula as obtained by the simulations, in particular the field is strongly compressed near the termination 
shock but with a lower than expected volume average. Ihe jet-torus structure is found to be clearly visible in high-resolution gamma- 
ray synthetic maps too. We also present a preliminary exploration of time variability in the X and gamma-ray bands. We find variations 
with time-scales of about 2 years in both bands. The variability observed originates from the strongly time-dependent MHD motions 
^^ ' inside the nebula. 



> 

m 



o 

OO 
O 






Key words, radiation mechanisms: non-thermal - Magnetohydrodynamics (MHD) - relativity - pulsars: general - ISM: supernova 
remnants - ISM: individual objects: Crab Nebula 



1. Introduction nism is, at most frequencies, synchrotron radiation by the rel- 

\^ ativistic particles gyrating in the nebular magnetic field. The 

_Cd Pulsar wind nebulae (PWNe, or plerions) are a class of super- synchrotron spectrum is cut off at some maximum frequency 

nova remnants (SNR), that originates from the interaction be- (typically in the soft gamma-rays) determined by the magnetic 

tween the ultra-relativistic wind blown by a pulsar (PW) and the field strength and the maximum energy that particles can attain 

surrounding supernova ejecta. The electromagnetic torques act- in the nebula. At higher photon energies the dominant emission 

ing on a fast spinning, highly magnetized neutron star convert mechanism is, most Ukely, inverse Compton (IC) scattering by 

its rotational energy into the acceleration of a cold magnetized the same particles interacting with different target photon fields. 

wind, expanding at relativistic velocity. The wind velocity needs For recent reviews of theoretical as well a s observational as - 

to be reduced in order to match the boundary condition of non pects of PWNe seeTOaensler & Slanel(l2006h : iKirk et al.1 (l2007h : 

relativistic expansion of the confining supernova remnant. This IBucciantinil (120061) . 

happens at a termination shock (TS), where the plasma is slowed rru * »j i . ■ .u . j tnwrt.^ 
, *^^ ,, , , • V ,.• ,-^ . . . , , The most recent developments in the study of PWNe are re- 
down and heated, the magnetic held IS amplihed and the nebular 1 t J » u- u u *• e ^u i^- ^ r\ ^u 

. , , , , , T n . . lated to high energy observations of these ob ects. On the one 

particle spectrum is thought to be produced. If conhnement is in- ,,.,,.. .• e r^u i j / » u u 

^ , ^ ^ ?.,,,•■ ,• . , , . hand, the latest generation of Cherenkov detectors has shown 

deed efncient, a non-negligible fraction of the energy lost by the mir^t . u .u i.- i, » • • 

, . , . . . ., , , • • , , , • , , , , PWNe to be among the highest energy emission sources in 

pulsar, which is invisible as long as it is locked in the cold rela- ^j^^ ^^^ ^ lAharonianlfeUOTl lGallant]|2007h . On the other 



tivistic Wind, may become detectable in the form of non-thermal . j j . i j .• i ■ ^ -^ c ■ j 

... ■ , , , .... .... , , hand, detailed spatial mapping at X-ray frequencies, as made 

radiation emitted by the relativistic particles in the nebula. ui u /-^u j u u ^u ^ ^u i- 

■' ^ possible by Chandra, has shown that the peculiar axisymmetric 

From the observational point of view, PWNe are character- morphology, known as jet-torus structure, initially seen in the 

ized by the emission of a very broad band spectrum of non- prototype of the c lass, the Crab Nebula ( Weisskopf et al. 200^ 

thermal radiation, typically extending from the radio to the X- iHester et al.ll2002l) . is also pre sent in a number of other objects: 

ray and even gamma-ray band. The primary emission mecha- the SNR associated to Vela dHelfand et al.l 120011: iPavlov et al.l 
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l2003h and to PSR B 1509-58 (iGaensleret ai]l2q02h : GO.9+0.1 
|Gaensleretal."2001'); G54. 1+0.3 dLu et alj2002l) : G130.7 + 3.1 
(Slane et al. 2004). Relativistic motions and time-varying small- 
scale features (wisps aroun d the TS and filaments in the torus) 
have also been reported (JWeisskopf et al.1 120001; iHesteret al] 
I2OO2I; iPavlov et al.1 l2003h . This paper intends to deal with all 
these aspects, presenting for the first time a study of the high- 
est energy emission from PWNe in the context of the jet-torus 
scenario, including an investigation of time variability. 

The theoretical interpretation of the jet-torus morphology 
observed in the optical and X-ray bands is based on the idea that 
the pulsar wind energy flux is anisotropic, depe nding on latitude 
above the pulsar rotati onal equator ( Bogovalov & KhangouliaiJ 
I2OO2I; lLvubarskvir2002l) . An energy flux that is maximum in the 
equatorial plane and decreasing towards the polar axis causes 
the ter mination shock to be oblate, closer to the pulsM" at th e 
poles dKomissarov & Lvub arskvl 120031: lDelZannaetalJl2004l) . 
Understanding the formation of a torus in the equatorial plane 
and of jets along the polar axis is straightforward in this scenario, 
where the post-shock flows first converge toward the equator and 
are then diverted along the symmetry axis by magnetic hoop- 
stresses. Recently, axisymmetric simulations in the relativistic 
MHD regime (RMHD) have been able to fully confirm t his view 
dKomissarov & Lvubarskvll2004t IPefZanna et al.ll2004b . 

The validity of this dynamical interpretation of the jet-torus 
structure was strengthened by the comparison between X-ray ob- 
servations of the Crab Nebula and non-thermal emission maps 
built on top of the flow struct ure resulting from RMHD simu- 
lations dDel Zanna et al.ll2006l) . In order to do this, information 
on the particle spectrum as a function of position in the nebula 
is needed. This kind of information is not available in a MHD 
approach, but Del Zanna et al. (2006) showed how the evolu- 
tion of the particle spectrum in the nebula, including synchrotron 
and adiabatic losses along streamlines, can be easily accounted 
for at least in an approximate way within a conservative MHD 
scheme, once an equation for the evolution of the maximum par- 
ticle energy is added to the code. The general method proposed 
in that work allows to compute synchrotron radiation, polariza- 
tion and spectral index maps based on MHD simulations. The 
resulting synchrotron surface brightness maps, when compared 
with the data, showed that the axisymmetric MHD scenario is 
not only able to account for the general jet-torus morphology, but 
also to explain a number of finer scale details. Optical and X-ray 
spectral index maps showed spectral hardening in the torus and 
softenin g toward the PWN border also comparable to t he obser- 
vations (iVeron-cettv & Woltieij[l993l:lMori et alj|2004 . Finally, 
the mildly relativistic nebular flow, with velocities in the range of 
those observed along the jets and around the TS (v ~ 0.5 - 0.8c), 
gives rise to features closely resembling the rings and the bright 
knot observed in the Crab Nebula, once Doppler boosting is 
properly taken into account. 

The combination of RMHD simulations and the diagnostic 
tools for synthetic synchrotron emission has thus proved to be 
a very powerful investigation tec hnique. Unknown paramet ers 
such as the wind magnetization cr (iKennel & Coronitilll984ah or 
the magnetic field angular distribution can be inferred by com- 
paring tJieJlieoredx^jiljrec^^ with observations. In particu- 
lar, in iDel Zanna et al.1 (|2006|) flie Crab Nebula optical and X- 
ray morphology were found to be best matched by assuming 
a latitude-averaged magnetization of (Tes = 0.02 and a nar- 
row striped wind region of low toroidal magnetic field along 
the equator (run A). The relativistic electrons produced at the TS 
were assumed t o be a single power law g' tza+i) ^jtjj q, - 0.6 as 
in the model bv lKennel & Coronitil(ll984bh . No attempt at mod- 



eling the radio emission was made. The resulting synchrotron 
spectrum was in excess of the data in the X-ray band (by a factor 
two) and show an unexpected flattening beyond a; 3 x 10'^ Hz. 

The aim of the present work is to further investigate the rea- 
sons for this difficulty in reproducing the observed high energy 
synchrotron spectrum. In particular we consider whether this 
could be solved through a different choice of the particle distri- 
bution function. Contrary to the magnetic field strength, the lat- 
ter does not affect the dynamics (provided a spectral index higher 
than 0.5) and can thus be changed with no consequences on the 
nebular morphology. However, while the synchrotron emission 
depends only on the combination of the magnetic field and parti- 
cle distribution, the IC emission allows to disentangle the contri- 
bution of the two. In order to exploit this fact in the present paper 
we extend the set of diagnostic tools for non-thermal radiation to 
the IC scattering process, and we apply the method to our simu- 
lations (runA) that best matched the Crab Nebula, extending the 
spectral coverage to the gamma-rays. 

As we mentioned the other recent development in the 
study of PWNe concerns the detection of very high en- 
ergy gamma-ray emission from these objects. TeV emission 
has been detected from the Crab Nebula (Aharoni anet al 
'2001 I2006dl: [Albert et al.ll2008l): MSH 15-52 (Aharonianet al 



120081): 
lerai] 



20053); Vela ( Aharonian s Tm l2006a[ Enomoto et al. 200^ 



Kookaburra complex (Aharonia n et all l2006b); HESS J1825- 
137 (Aharonian etal. 2 006c) ; the composite SNR G 0.9+0.1 
(Aharonian et al. 2005a! ); die two candidates H ESS J 1357 - 645 
and HESS J1809- 193 (lAharonian et al.ll2007l) . In this context, 
TeV emission can in principle result from two different pro- 
cesses: either inverse Compton emission involving relativistic 
electrons that up scatter lower energy photons, or n'-'' decay in- 
volving the presence of relativistic protons that produce pions by 
nuclear collisions. The presence of relativistic protons in PWNe 
is also suggested by theoretical reasons. In fact, particle accel- 
eration at transverse relativistic shock waves, such as the pul- 
sar wind TS, is difficult to explain through standard acceleration 
processes. The most successful model so far is based on resonant 
absorption by electrons and positrons of the relativistic cyclotron 
radiation produced by ions, that are predicted to be present in the 
pulsar wind. The model also requires the ions to be energetically 
dominant in the wind, though very few by number. 

The most likely mechanism at the origin of this emission is 
inverse Compton scattering by the same particles that produce 
the nebular synchrotron spectrum. These particles can up scatter 
different targets: CMB photons, synchrotron emitted photons, 
starlight and possibly FIR photons due to reprocessing of the 
nebular radiation by dust within the SNR. Studying IC emission 
is important for two different reasons. First of all, as anticipated, 
modeling synchrotron and IC radiation at the same time allows 
to disentangle the information on the magnetic field strength and 
particle number density which is always combined w hen consid- 
ering synchrotron emission alone (Gould 1965; De Jager et al.l 
1996). Moreover, a detailed modeling of the IC component is at 
present the only way of investigating whether the high energy 
data leave room fo r an extra contribution of had ronic origin (e.g. 
lAmato et al.ll2003l: iBednarek & Bartosikll2003l) . Clarifying this 
point is also important in view of the quest for galactic sources 
of cosmic rays at energies around the knee. 

This paper presents a general method that could be applied to 
the entire class of PWNe and more generally to all non-thermal 
sources. In fact, it consists in evolving in time and space the max- 
imum energy of emitting particles, together with the other MHD 
dynamical variables, while the appropriate simulation set up 
(e.g. initial and boundary conditions, resolution) and the shape 



D. Volpi et al.: Non-thermal emission from relativistic MHD simulations of PWNe: from synchrotron to inverse Compton 



of the distribution function at injection sites are arbitrary and 
can be both chosen as the most appropriate for the object un- 
der investigation. Here, however, our purpose is not to optimize 
the dynamical free parameters to reproduce specific nebulae, but 
to apply our emission model to the runA simulation and com- 
pare the results to the Crab Nebula data. The Crab Nebula is, 
in fact, the brightest PWN also at very high energies, and is re- 
garded as a standard candle for gamma-ray observations. This 
makes it a natural target for new instruments and a wealth of 
data is already available from ground and space instruments: 
EGRET (Nolan etal. 1993), COMPTEL (Kuiperetal. 2001), 
HEGRA (Aharonian et al. 2004), HESS ( Masterso n et al.1120051: 
lAharonia n et al. 2006d), MAGIC (Albert et al. 2008]) . This year 
GLAST will be launched and so new data will be available from 
20MeVto300GeV. 

The existing gamma-ray observations are here used to con- 
strain the parameters of our model. The synthetic emission 
up to TeV energies is calculated for the first time based on 
time-dependent 2-D numerical simulations. In particular we 
assume the configuration coiTesponding to runA, as reported 
in [Del Zanna et al. (2006), and following Atovan & Aharonian 
d 1 9961) we consider a more general distribution function of emit- 
ting electrons consisting of two families: one corresponding to 
radio emitting electrons, that can be interpreted as primordial 
population of particles and the other coiTesponding to particles 
continuously accelerated at the termination shock, responsible 
for the synchrotron spectrum up to the gamma band. The IC 
emissivity is calculated (in the optically thin regime, appropri- 
ate for these objects) using the general Klein-Nishina differen- 
tial cross-section. The three different photon targets recognized 
as the most important for the Crab Nebula (Aharonian et al.l 
l2006dl) are considered: the nebular synchrotron photons, those 
respo nsible for the far-infrared thermal excess (iMarsden et al.l 
1 19841) and the cosmic microwave background. Previous work 
applied similar emission reci pes only to stationary and radi- 
ally symmetric models (e.g. lGouIdlll979.), and notic eably to 
the Ke nnel & Coroniti (KC: iKennel & Coronitilll984ah RMH D 



model dPe Jager & Harding|ll992HAtovan & Aharonianll996l) . 
The paper is structured as follows. In Sect.|2]our non-thermal 
emission model is presented. The integrated spectra as arising 
from both the KC model (used here to test the emission recipes) 
and our simulations are shown in Sect.[3]and compared with ob- 
servations of the Crab Nebula. In Sect. |4] we produce surface 
brightness maps at different energies in the gamma-rays to com- 
pare with existing and future images. Finally, time variability is 
studied in Sect. |5] linking the synthetic X and gamma-ray emis- 
sion. 



2. The non-thermal emission model 

In this section the formulae used to compute synchrotron and 
IC emission are explained in detail. As mentioned in Sect. [1] 
the evolution of the distribution func tion and the syn c hrotro n 
emission recipes are adapted from iDel Zanna et alj (l2006h . 
while the choice of the distribution function at the TS follows 
lAto van & Aharonian (1996). Some approximations are neces- 
sary to avoid using formulae in which some terms are unknown: 
for instance it is not possible, within the present scheme, to as- 
sociate the local pressure value to the coiTesponding one (along 
streamlines) at the termination shock. This would require evolv- 
ing in space and time the entire distribution function, a task 
which is beyond the goal of this paper. As we will see a few 
additional approximations are used with the aim of reducing the 



computational costs, after checking that they do not change sig- 
nificantly the final results. 

2.1. Relativistic electrons 

In the Crab Nebula spectrum two main breaks appear: one in the 
IR and the other in the UV. In order to account for both breaks 
we consider the presence of two distinct populations of relativis- 
tic particles (electrons and positrons). In analogy with previous 
MHD models we assume the break in the UV to be due to syn- 
chrotron cooling of a population continuously accelerated at the 
TS and responsible for the o ptical and higher frequency radiation 
dKennel & Coronitlll984bt) . We then describe the radio emitting 
electrons as a different popul ation, which has sometime s been 
interpreted as primordial (e.g. 'Atoyan & Aharonian' 1996^, born 
at the outburst of the supernova. The distribution function of b oth 
populations is shaped following lAtovan & Aharonianid 19961) . 

Assuming isotropy, the radio emitting particle distribution 
function per unit solid angle is taken as 



Me) = 



4^' 



(2ff,-+I) 



exp(-e/e;). 



(1) 



assumed to be homogeneous in the PWN and constant in time. 
Here e is the normalized energy of the particle (its Lorentz fac- 
tor), Al is the normalization constant, Oi- is the radio spectral 
index, e* is the radio cut-off energy, which is the energy cor- 
responding to a synchrotron cooling time comparable to the age 
of the nebula. 

The second population is composed by the relativistic wind 
particles continuously accelerated at the termination shock and 
injected downstream into the nebula, with the following distri- 
bution function 



/o(eo) = /(e + eo)"<'""^" exp(-6o/e:). 
47r 



(2) 



This has to be evolved along streamlines in the PWN and may 
change with time. The variables with the index are calculated 
at the TS. As for the previous case Aq is a constant to be deter- 
mined based on the local number density or pressure, a^, is the 
wind spectral index and e^ is the cut-off energy (corresponding 
to the maximum synchrotron frequency in the gamma-rays). The 
additional constant e represents the minimum energy of the wind 
population, and will be chosen to match smoothly the infrared- 
optical synchrotron emission of the two families of electrons. 

As far as the evolution of the wind population is concerned, 
following iDel Zanna et al.l d2006l) . adiabatic and synchrotron 
losses are taken into account by defining the evolved distribu- 
tion as 



4/3 



/^(,) = (^j (^)V„(eo), 



where p is the rest mass density, thus: 



/w(e)=^P 



(?)< 



e + eo) 



-(2ff„ + I) 



exp (-eo/e;). 



(3) 



(4) 



Here we have used the conservation of particles' number along 
streamlines and assumed adiabaticity for the pressure (p oc p'*^-'). 
A() = A^po, and we take into account synchrotron losses through 
the quantity eoo, which is the maximum possible particle energy 
at a position along the streamline: 



eo 



1 - e/e« 



(5) 
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The variables p, p and e^o are in general functions of space and 
time and they will be here provided by our code. Notice that 
adiabatic losses are taken into account only through the pres- 
sure term in Eq. (HJi, while we are forced to neglect them else- 
where (terms (p/po)'^' * are approximated to 1), as explained in 
iDel Zanna eFd] (l2006l) . 



2.3. Inverse Compton emission 

The Inverse Compton spectral power emitted by a single ultra- 
relativistic particle with normalized energy e may be written as 



rf(y,E) 



""/(^b 



(v,, V, e) «v(v,) dvt, 



(12) 



2.2. Synclirotron emission 

The synchrotron spectral power emitted by a single ultra- 
relativistic particle is 



Pl™{v, e) = 2cTjc-^e^ 6(v - y„), 
on 



(6) 



where o-j is the Thomson cross-section, B± is the local mag- 
netic field component normal to the particle's velocity and 
V is the observed frequency. Notice that for simplicity we 
have also considered the monochromatic approximation, where 
Vm - 0.29vc is the location of the maximum emission (see e.g. 
iRvbicki & LightmanI 19791) and the critical frequency is given by 



3e 
47Tmc 



B.e' 



The emission coefficient is, in the general case 



jfyv)^^'pf\v,e)f{e)Ae, 



(7) 



(8) 



whatever the distribution function /(e). Thanks to the monochro- 
matic approximation we can readily perform the above integra- 
tion and find 



Jv (y) = 2o-tc— :; -s-f{e{v)\. 



2v Stt" 



(9) 



where 



where tiyivi) is the number density of target photons per 
unit frequency v^. The general form of the differential cross- 
section for IC scattering per unit frequency, incl uding both th e 
Thomson and Klein-Nishina regimes, is given bv lJonesI (Il968l) : 
iBlumenthal & Gouldl(ll970h : 



dcr\ 3 cTj 



2^1n^ + (l -^)|1 +2q+-^^^^'' 



2 1 +r^ 



with 



q 



hv/mc 



r(e — hv/mc^) 



and r = Aehvilmc . 



(13) 



(14) 



Due to the kinematics of the scattering, q' is a parameter in the 
range < q' < 1 and F determines the energy domain (the 
Thomson limit corresponds to F <K 1). 

As in the previous case, the emission coefficient is given by 



7f(v) 



-S" 



n^(y,e)/(e)de. 



(15) 



for any distribution function of scattering particles /(e). 
However, in the most general IC scattering case we cannot 
reduce the above formula to an analytical expression and a 
double integral has t o be computed numerically. Following 
IBlumenthal & Gouldl (Il970l) . we replace the integration over e 
by an integration over q. The relation between these two vari- 
ables is 



^•^^^ = ~I3 



1 hv 

2 mc^ 



1 + 



\ + sq 
sq 



1/2 



hv hv, 
^ ,, (16) 



e(y) = |0.29-^B^ 
Amnc 



-1/2 



.1/2 



(10) 



It is easy to verify that, in the power law regions of both the dis- 
tribution functions described in the previous section, the emis- 
sion coefficient behaves as yV oc v^" , as expected. 

Due to relativistic beaming, the magnetic field component in 
the above expressions can be written as B^ = |n x B|, where n is 
the line of sight direction. To reduce computational costs, here 
we will consider an isotropically distributed field (discrepancies 
in the integrated fluxes are less than 3% at 1 keV), for which 
B^ - ■\/2/3B/'y, B is the toroidal magnetic field as measured 
in the laboratory frame and y is the Lorentz factor of the fluid 
element. The synchrotron spectrum, expressed in terms of the 
total luminosity per unit frequency, is 



^^(v) 



/4.;t™(r, 



v)dV 



(11) 



where the integral is over the nebular volume. 

Synthetic surface brightness maps are derived instead by 
integrating the emissivity along the line of sight alone. The 
Doppler boost ing towards the observ er is taken in account as 
in the paper bv lDel Zanna et all (l2006h . 



where the dimensionless quantity s is small in the Thomson Umit 
and large in the extreme Klein-Nishina regime. 

The target photon density in Eq. (fT2l i may have different ori- 
gins. Here we study three cases, recogni zed to be relevant for 
the Crab Nebula and for PWNe in general (lAtovan & AharonianI 
1996): IC scattering of reprocessed synchrotron photons (IC- 
SYN), IC scattering of the far-infrared thermal radiation by lo- 
cal dust (IC-FIR), and IC scattering of cosmic microwave back- 
ground photons (IC-CMB). 

In the self-synchrotron emission case, the photon den- 
sity at any position r can be determined from the radiative 
transfer equation in the optically thin re gime (G ouldl 119791 : 
iDe Jager & Hardin d[T992t lAtovan & Aharonian. 19961) : 



«'"-'™(r. 



chv, j 



;SYN 



(r',Vr) 



|r' 



dV\ 



(17) 



where the spatial integral is over the whole nebula. However, in 
order to reduce the computational costs, after testing that differ- 
ences are negligible, we use here the approximation of spherical 
symmetry and of a homogeneous synchrotron emission coeffi- 
cient. We thus replace the above formula with 



«r'™(^n) 



1 ^^^^(v,) 
chv, AnR^ 



U{r/R), 



(18) 



D. Volpi et al.: Non-thermal emission from relativistic MHD simulations of PWNe: from synchrotron to inverse Compton 



where the synchrotron spectral luminosity is given in Eq. ( fTTT i. 
R is the radius of the nebula and 



2 Jo X \x-y\ 



with X - r/R and y - r JR. U( x) is a non-dimensional quan tity 
decreasing with x from 3 to 1.5 JAtovan & Aharonianll 19961) . 

When the emission recipes are applied to the results of our 
simulations, the synchrotron photon density in Eq. ( fTSl l becomes 
a function of time. Variations of the synchrotron emitted flux 
are up to 30% on time-scales of order a few years, as we will 
discuss later on. Since this time-scale is not much larger than 
the light crossing time of the inner nebula, the temporal lag be- 
tween the emission of a photon and its comptonization should be 
taken into account. This could properly be done only by using 
more complex, time dependent radiation transfer equations, that 
should be solved simultaneously with the MHD dynamics. Such 
an analysis is beyond the goal of the present paper, where we will 
make the approximation of instantaneous comptonization. This 
implies that our results on gamma-ray variability at frequencies 
involving a time-varying photon target density should be taken 
with caution. 

In the case of the thermal FIR emission, the emitting dust 
is expected to be located in a torus created by the sup ernova 
progenitor and in the optical filaments JGreen et al.ll2004t) . Since 
the dust volume and location are uncertain, we choose the crude 
approximation: 



nf-^'\v,)^ ^ ^^ ^""'^ 



c hvt 4nR^ 



(20) 



with L™ corresponding to a uniformly emitting region of black- 
body radiation. 

Finally, the CMB photon density is provided by the Plank 
distribution 



, IC-CMB 



(Vr) 



Sn 



c^ expihvr/ksT) - 1 



where k^ is the Boltzmann constant and T -2.7K. 

Once the IC emission coefficient has been calculated as a 
function of space, integration over the PWN volume gives the 
total spectral luminosity as in Eq. ( fTTT i 



l:"( 



-/ 



47r;f(r,v)dy, 



(22) 



where scattering with different targets will be computed sepa- 
rately to investigate the relative contribution to the overall IC 
emission. 

The surface brightness maps with the relativistic corrections 
will be computed as for synchrotron emission. 

3. Integrated spectra: fitting the Crab Nebula data 

In order to understand the different contributions of the various 
physical processes to the overall emission of the Crab Nebula, 
simulated synchrotron and IC integrated emission spectra are ob- 
tained and plotted against the data from the radio to the gamma- 
ray band. 

From this comparison, we will derive the unknown parame- 
ters in the distribution functions defined in the previous section, 
within two different dyna mical models: the steady-sta te, spheri- 
cally symmetric model bv lKennel & Coronitil(ll984al) . used as a 
test, and our axisymmetric 2-D simulations. 



Table 1. Values of the parameters for the spherical model (KC) 
and for two different set up of our axisymmetric RMHD simula- 
tions, runA (1) and runA (2). Ar is in unit cm""* and A„ is in unit 



(19) erg- 



Parameters 


KC 


RunA (1) 


RunA (2) 


a. 


0.26 


0.26 


0.26 


OTw 


0.70 


0.70 


0.85 


Ar 


2.90 X 10-" 


1.35 X 10-= 


1.35 X 10-= 


A^ 


3.2 X 10** 


2.6 X 10' 


3.2 X 10" 


< 


2.0 X 10= 


3.0 X 10= 


3.0 X 10= 


< 


1.4 X 10'° 


oo 


oo 


6 


4.3 X 10= 


5.0 X 10= 


8.0 X 10= 



3. 1 . The spherical model 

The emission recipes of Sect. |2] are here applied to Kennel and 
Coroniti's spherical model, c alculated for a win d magnetization 
parameter cr - 0.005 (as did lAtovan & Aharon ian 1996) and a 
nebular radius R - 2.0 pc. Integrated spectra in Fig. [1] are com- 
pa red with observations of the C rab Nebula and with the results 
by lAtovan & AharonianI (Il996l) . used here as a benchmark for 
our model. 

The values of the parameters correspond to best fit es- 
timates (see Table [TJ. The primordial and wind spectral in- 
dices are equal to the observed values in the radio and optica l 
bands dBaars & HartsuiikeiJll972l:IVeron-cettv & WoltieiJll993h . 
respectively. Ai- and Aw are found by fitting the emission in the 
radio and X-ray band; e* and e are chosen so as to connect the 
radio and optical parts of the spectrum smoothly; e^ is such as 
to reproduce the observations in the gamma-ray band. We re- 
call here that Aq - A^/Po refers only to the non-thermal tail 
of the particle energy distribution. This accounts for some un- 
known fraction f of the th ermal pressure at the termination shock 
JKennel & Coronitilll984bl) . with ^ determined from: 



/TIN mc^ r r*^™ 

(21) ^^0 = ^- d^ /o(eo)eodeo 0<^<1, (23) 



where Q. is the solid angle, e™" - e, e™"" = e^ and fo(eo) con- 
tains Aq and hence pQ. We find ^ * 70%. 

The values of o ur parameters are sli g htly d ifferent from those 
in the paper by Atova n & AharonianI (1 1 996*). from which the 
form of our distribution function is taken. The main difference 
is in the cut-off energy e^ and is due to our approximated for- 
mula (|5]l, where adiabatic losses are neglected (these are con- 
sidered instead by means of the factor p in Eq. (|4]l). However 
the difference involves only the synchrotron spectrum at gamma- 
ray energies, with no consequences for the IC radiation flux, to 
which electrons near the cut-off contribute negligibly (see also 
Fig.|6]and the discussion in Sect. |5]l. 

In the left panel of Fig. [T]we plot the total evolved distri- 
bution function, given by Eqs. ([TJ and (|4|, for different values 
of eco. The particle distribution function is normalized through 
a value of the local pressure corresponding to the volume aver- 
age in the simulation. Since the relic population is assumed to be 
spatially homogeneous, the radial dependence is only included 
in the wind particle distribution function thro ugh the local values 
of p(r) and eoo{r), both provided by the model lKennel & Coronitil 
(Il984allbl) . The wind population displays the expected behavior: 
both the local value of the thermal pressure, to which the peak 
of the distribution is proportional, and the cut-off energy, pro- 
viding the exponential decay, are decreasing functions of r. Our 
results are basically coincident with those shown in Fig. 3 of 
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Fig. 1. Left panel: Kennel&Coroniti total electron energy distribution function for different values of the maximum particle energy 
eoo as specified in the figure. In the plot we use E - emc^ and F - (f^ + /„) ■ Anl{mc^) and the particle distribution function is 
normalized to an average value for the pressure of 2.1 x lO^^erg cm"^. The dashed lines are the primordial and the wind component 
at eoo — » °°- Right panel: Kennel&Coroniti spectra are plotted against the data. The single contributions from different incident 
photon targets, thermal FIR emission and total luminosity are reported with different colors and line-styles. Dashed lines: primordial 
and wind synchrotron contributions and IC-SYN emission; dash-dotted line: thermal FIR and IC-FIR ones; dash-dot-dotted line: 
IC-CMB one; solid line: total luminosity. References to the observations reported here are given in the main text. 



JAtoyan & AharonianI d 1 9961) in spite of a different visualization 
(see the respective captions). 

The overall synthetic spectrum is then calculated using 
the above distribution function with the dynamics taken from 
Kennel&Coroniti's model. In the right panel of Fig.[T]the spec- 
tral luminosity is plotted and the individual synchrotron contri- 
butions of the two populations are also shown. The synchrotron 
spectrum from radio to IR frequencies is due to the primor- 
dial particles, homogeneously distributed all over the nebula. 
At higher frequencies the emission is due to the wind parti- 
cles, which suffer from adiabatic and synchrotron cooling. The 
FIR thermal radiation is obtained from the black-body for- 
mula with a temperature of the emitting du st of T = 46 K 
(JStrom & GreidanusllI992l; iGreen et alj[2004 . As far as IC ra- 
diation is concerned, we find that the primordial electrons con- 
tribute to IC-CMB radiation up to 1.6 x 10^^ Hz, to IC-FIR up to 
3.2 X 10^3 Hz and to IC -SYN up to 7.9 x lO^^ Hz (in the paper 
by JAtovan & Aharoni an 1996, the frequency value is larger by 
almost one order of magnitude) mainly in the Thomson regime, 
ehvt <c mc^, while at higher frequencies the emission is due to 
scattering by the wind particles. 

References for the reported d a ta are as follows: radio 
data are from iBaars & Hartsuiikerl (|I972|); th e mm data are 
from iMezger et all (Il986h and Ifiandiera et alj (|2002|) ; the in- 
frared points are from I RAS in the far to mid-infrared 
dStrom & GreidanusI 1 19921) and from ISO in the mid to 
near infrared range JDouvion et al.l 1200 lb ; optical is from 
IVeron-cettv & Woltieij (Il993l) and UV from JHennessv et~an 
(|I992|) . Points in the range betwee n soft X and gamma-rays 
are taken from 'Kuipe r et al.l (|200l'), who compiled data from 
BeppoSAX, COMPTEL and EGRET. In the TeV band we 
plot the data from M AGIC jAlbert et alj 120081) and HEGRA 
JAharonian et alj|2004 . 

Shapes and values of the single curves, especially IC- 
SYN and IC-CMB , are consistent with those of Figs. 9 of 
JAtovan & Aharonia n ( 1996), so we conclude that in our model 
the computation of inverse Compton emission seems reliable. 
We recall that the main simplifications introduced here are in the 



wind distribution function, Eq. (|5]l, and the homogeneous emis- 
sivity in the IC-SYN target photon density, Eq. ( fTSl ). 

The computed synchrotron emission is in agreement with 
the data, whereas the IC luminosity is lower at all gamma- 
ray frequ encies. In both the present 1 -D calculation and in the 
model by JAtovan & AharonianI d 1 9961) the gamma-ray EGRET 
and TeV data can be simultaneously fitted only by invoking addi- 
tional contributions to the emission namely Bremsstrahlung and 
hadronic processes. In addition one may tune the magnetic field 
strength, however this has only been done in an ad hoc manner 
(JAtovan & AharonianI! 1 9961; JAharonian et alj|2004 rather than 
changing the dynamics self-consistently. 

The conclusion that the IC alone is insufficient to explain the 
gamma-ray fluxes from the Crab Nebula seems to be the com- 
mon result of all 1-D models. However, 1-D stationary models 
are not able to account for a number of observed nebular prop- 
erties, in particular the jet-torus morphology seen in the X-rays. 
Due to the oversimplified dynamics intrinsic to spherical mod- 
els, the conclusions about the emission processes are also ques- 
tionable. On the other hand, 2-D time-dependent RMHD models 
have proved to reproduce the inner structure of the Crab Nebula 
even in some of its finest details. In the following sections, we 
present a study of the synchrotron emission and, for the first 
time, of the IC gamma-ray radiation from PWNe in the context 
of the axisymmetric scenario. 

3.2. Results from numerical simulations 

3.2.1. General discussion 

For the study of the non-thermal emission in the context of a 
2-D model we focus here on the simulation labelled as runA in 
the paper byl Del Zanna et all (I2006.) . The parameters that deter- 
mine the shape and the strength of the magnetic field are: the 
anisotropy parameter, a = 0.1, the width of the striped wind re- 
gion, b = 10 (corresponding to » 10° around the equator), and 
the wind magnetization, cr - 0.025. The particles' maximum 
energy at injection, too, is here taken as 10'", corresponding to 
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Fig. 2. Left panel: runA ( 1 ) total electron energy distribution function versus particle energy for different values of maximum particle 
energy foo as specified in the figure. In the plot we use E - emc^ and F - (f^ + /„) ■ Anjimc^) and the particle distribution function 
is normalized to an averaged pressure value of 2.1 x 10 ''erg cm"^. The dashed lines are the primordial and the wind component for 
foo — » oo. Right panel: runA (1) spectra are plotted against the data. The single contributions from different incident photon targets, 
thermal FIR emission and total luminosity are reported with different colors and line-styles. Dashed lines: primordial and wind 
synchrotron contributions and IC-SYN emission; dash-dotted line: thermal FIR and IC-FIR ones; dash-dot-dotted line: IC-CMB 
one; solid line: total luminosity. References to the observations reported here are given in the main text. 
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Fig. 3. The same as Fig.|2]for runA (2). 



5x10 ^^ eV (see e.g. JDe Jager & Hardindl992l;lAharonian et al.l 
12004 . 

We now discuss the nebular emission spectra resulting from 
our simulation. Two different distribution functions are consid- 
ered for the wind particles: one with the wind spectral index 
ttw = 0.7, runA (1), and the other with a„ - 0.85, runA (2). 
The first one is the counterpart of the spherical model computed 
in the previou s section, and coiTesponds to the best fit optical 
spectral index dVeron-cettv & WoltieilI993h . However, with this 
value of cKw we are far from fitting the synchrotron spectrum 
from X to gamma-ray frequencies and the IC emission exceeds 
the data, as it is clear from the right panel of Fig. |2] A steeper 
injection spectrum is then adopted in order to discuss how the 
combined constraints from synthetic synchrotron and IC emis- 
sion can help determining the physical parameters of PWNe. 
This second spectrum, in fact, allows one to reproduce the high 
energy synchrotron spectrum (see Fig. [3]) but the TeV emission 
is slightly closer, though still above the data. In both runs the 



spectrum around lOGeV seems approximately to fi t the EGRET 
data, that, we recall, are affected by large error bars (' Kuiper et al.l 
1200 ll) . Hopefully GLAST observations will allow to constrain 
better the models in this crucial range of frequencies. At the mo- 
ment, however, there is no room for additional contributions to 
the gamma-ray emission, such as Bremsstrahlung and hadronic 
processes. Note that these results are opposite to what found in 
1-D (see Sect. 13. Il l, where the model underpredicts the IC peak 
emission of a factor =» 5. 

Before discussing the spectra in detail, let us first describe 
the particles distribution function. The total (primordial plus 
wind) particle distribution functions for the two cases, obtained 
as in Sect. 12.11 are shown in the left panels of Figs.|2]and[3] The 
values adopted for the maximum electron energy too and for the 
thermal pressure are the same as those used in Fig. [T] However 
one should remember that these quantities are now functions of 
both the spherical radius r and the polar angle 6 (see Fig. |4|i. 
The remaining parameters are chosen such as to obtain with our 
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Fig. 4. Left panel: color coded map of the thermal pressure in units of erg cm"^. Right panel: color coded map of the maximum 
particle energy; superimposed contours correspond to values of 10^ (solid line), 10^ (dotted line), 10^ (dashed line) and 10'" (dot- 
dashed line). Along the axes, distance from the central pulsar is reported in ly. 



simulated spectra the best fit possible of the observations (see 
TableQ. 

By comparing the left panels of Figs. |2] and [3] with the cor- 
responding panel of Fig. [T] one can see that in our simulations 
a larger number of electrons is present at all energies. However, 
in runA (2) the distribution function falls down rapidly at high 
energies reaching values that are close to those found within the 
spherical model. 

The values of the parameters adopted are reported in Table[T] 
Notice that for the assumed injection energy, e^ is never larger 
than 10'". Hence the synchrotron cut-off in the gamma-rays is 
coiTectly reproduced without the need of introducing an expo- 
nential decay of the particle distribution function. 

For both cases, runA (1) and runA (2), the emission spec- 
trum is calculated as described in Sect. |2] and displayed in the 
right panels of Figs. |2] and [3] As far as the different contribution 
of the two populations in different frequency bands is concerned, 
the primordial electrons contribute to the synchrotron emission 
up to IR frequencies, to IC-CMB emission up to 2.5 x 10^^ Hz, 
to IC-FIR up to 3.2 X 10^3 Hz and to IC-SYN up to 7.9 x lO^^ Hz 
mainly in the Thomson regime. The radiation at higher frequen- 
cies for all components is due instead to the wind electrons. In 
both runs, the data are well fitted in the radio (by construction) 
and in the optical bands. In this context the first spectral break, 
between the IR and the optical, appears to be intrinsic and due 
to the superposition of the two different populations of parti- 
cles. The second one, in the UV band, corresponds in our model 
to a strong synchrotron burn-off and is related to the decreas- 
ing volume occupied by particles of increasing energy. Beyond 
Si 2 X 10'^ Hz, we observe a spectral flattening in our curves, 
followed by further minor break s and counter-breaks at hig her 
frequencies. Here, we recall that Kargaltse v & Pavlov! ( 1200 8') re- 
port X-ray observations of PWNe with spectral indices less than 
0.5, which have no theoretical explanation but that in principle 
could be reproduced by simulations similar to runA (1). 



In the case of runA (1), at frequencies above x 10'* Hz the 
computed synchrotron emission is always far in excess of the 
data. The same happens at frequencies larger than » lO^'' Hz, 
where inverse Compton is the dominant emission mechanism. 
Notice however that the IC-SYN component is not much af- 
fected by the hard X and gamma-ray synchrotron excess, which 
contributes only beyond 1 TeV. In runA (2) we manage to fit the 
X-ray data with our steeper simulated spectra, but still the IC 
emission appears to exceed the data by a factor ~ 2. 

These results show that the difficulties we have in reproduc- 
ing the synchrotron spectrum with the commonly accepted value 
of a„ cannot be solved by simply changing the particle distribu- 
tion function, but are rather hinting to a problem related to the 
underlying dynamics, and in particular with the magnetic field 
distribution in the nebula. 

In our simulations the nebular magnetic field appears to be 
compressed in localized regions close to the termination shock, 
while rapidly decreasing outward, much faster than in spherical 
models. Its volume averaged value is Bmean = 1 x lO"'* G, which 
is 2 - 3 times less than the value inferred from previous models 
dKennel & Coronitill984al) . This means that in the outer regions, 
which occupy a large fraction of the volume and contribute pro- 
portionally to the integrated emission, our synthetic magnetic 
field is much lower than for spherical models. Hence, in order 
to fit the integrated spectra, we were forced to introduce a larger 
number of emitting electrons compared to the spherical model. 
Our electrons are also subject to less severe synchrotron losses 
on average, and consequently a steeper spectral index a„, as in 
runA (2), is necessary in order to reproduce the high energy syn- 
chrotron spectrum. The simultaneous calculation of the IC spec- 
trum, however, shows that the required large number of emitting 
electrons is not allowed by the gamma-ray data. Therefore, in 
order to reproduce the spectrum of the Crab Nebula, a larger 
and more distributed magnetic field would be needed, provided 
that it is able to reproduce the correct dynamics. Nevertheless 
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before drawing any conclusions on the physics of the emission 
processes at work, a further investigation in the parameter space 
(for both the dynamics and the particle distributions) should be 
first attempted, since here we have only considered the runA re- 
sults. Possibly the problem of the concentration of the nebular 
magnetic field is intrinsic to axisymmetric simulations and could 
be alleviated only by moving to 3-D. 

3.2.2. On the spectral slopes 

Let us now go back to the most surprising features observed in 
our spectra, namely the multiple changes of slope, which are 
in contrast with simple expectations from synchrotron cooling. 
However, the prediction of only one single steepening of the 
slope due to synchrotron burn-off is based on simplified 1-zone 
or 1-D nebular models, that assume a smooth flow and magnetic 
field distribution, so that particles experience continuous syn- 
chrotron losses. On the contrary, one of the main results of 2-D 
MHD models of PWNe has been the realization that the inter- 
nal flow dynamics is more complex than previous 1-D models 
suggested. 

This complexity reflects on the particle spectrum as a func- 
tion of position within the nebula. As anticipated this is com- 
puted as described in Sect. 12. 11 using the information on the ther- 
mal pressure p and on the maximum energy eoo provided by the 
simulations: color coded maps of p and eoo are shown in the left 
and right panel of Fig. ^respectively. 

The computed pressure tends to a constant with a spatial av- 
erage of 2.1 X 10^'erg cm"^ (which is behind our choice of the 
normalization constant in the plots showing the population dis- 
tribution functions in Figs. [T]|2] and |3j. Exceptions are observed 
along the jets, where the maximum is located, and in the equato- 
rial plane, where vortices create regions of lower pressure. 

As far as the local value of the maximum energy is con- 
cerned, the effects of synchrotron burn-off are highlighted by 
the four contours added on top of the color coded map of eco . 
These correspond to the same values for which particle distribu- 
tion functions are shown in the left panels of Figs.[Tl|2]and[3] 
The main difference between the results of our 2-D simulations 
and spherically symmetric models is immediately apparent from 
the contours in the right panel of Fig. [U the spatial domain oc- 
cupied by particles with energies up to a certain value is made of 
disconnected regions. While in 1-D the contours would be con- 
nected and encircle progressively smaller regions of space with 
increasing values of e^o, here the flow structure and magnetic 
field distribution is such that this is not the case anymore. In fact 
in our simulations fast flow channels in regions of relatively low 
magnetic field may take high energy particles far from the ter- 
mination shock without suffering important synchrotron losses. 
In addition to this, the spatial distribution of the magnetic field 
arising from our model, which is compressed b y the flux vortices 
aroun d the TS and the polar axis (see Del Z anna et al] |2004 
|2006|) . causes a non monotonic decline of the volume occupied 
by high energy particles. This is at the origin of the multiple 
changes of slope in the synthetic emission spectrum. In fact, it 
is interesting to notice that 2-D simulations performed with a 
lower value of the magnetization parameter give origin to much 
less structured flow patterns and spectra showing just one change 
of slope, analogous to those found from spherical models. At the 
same time, however, the jets disappear. On the other hand, based 
on the possibility we find here of having more than one spec- 
tral break due to synchrotron losses, one may speculate that a 
larger value of the magnetization parameter, sufficient to move 
the multi-slopes to IR frequencies, may allow to explain the en- 



tire synchrotron spectrum of the Crab Nebula with just one popu- 
lation of particles. However additional particle acceleration may 
take place in regions other than the termination shock, in partic- 
ular the termination spots of high speed flows in the equatorial 
plane and along the jets, and gives rise to the changes in slope. 

3.2.3. On the SYN/IC estimate for B„ean 

The ratio between synchrotron and IC observed luminosity, 
when these data are available, is often used to evaluate the aver- 
age magnetic field strength in PWNe (e.g. De Jager et al. 199g; 
IH. E. S. S. CoUaboration: A. Diannati-Atai et al.ll2007l and ref- 
erences therein). However such an estimate is based on the as- 
sumption that the magnetic field is uniformly distributed (or 
slowly varying as in the KC model). On the other hand the re- 
cent numerical results of PWNe have clearly shown (as also con- 
firmed in this paper, see the discussion in this section) that the 
distribution of the magnetic field and of the maximum energy of 
the emitting particles can change quite dramatically within the 
nebula. 

It is thus of interest, not just to observers but also to theorists, 
to understand if, in the case of 2-D PWNe models, the magnetic 
field inferred from the ratio SYN/IC provides a good approxima- 
tion of the true value, or a biased one, at different frequencies. 

In particular we compare the estimate of the magnetic field 
from the ratio SYN/IC, that are computed from the synthetic 
spectra, with the average value of the magnetic field Bmean equal 
to 1 X 10 '^G in our RMHD simulations, on which the same 
spectra are based. The obtained discrepancy is of order 10% at 
synchrotron frequencies below 10^^ Hz, while beyond 10^^ Hz 
the discrepancy becomes larger than a factor of 2. This is due 
to the fact that in the high frequency range the emission comes 
from wind particles confined to a small volume around the TS, 
where the magnetic field is highly structured and higher than the 
average value. 



4. Gamma-ray surface brightness maps 

In the present section, synthetic IC surface brightness maps are 
shown, for the first time, in the gamma-ray band. These are 
calculated for runA (2) by integrating along the line of sight 
the total (IC-SYN, IC-FIR and IC-CMB) emissivity, where the 
Doppler relativistic effects are included as explained in the paper 
bv lDel Zanna et all (I2006I) . 

In Fig. |5] images at photon energies of 4GeV, 250 GeV and 
1 TeV are computed. These are chosen to lie in the middle of the 
band which wiU be observed by GLAST (20MeV-300GeV), in 
the range measured by MAGIC (60 GeV-9 TeV), and at the cen- 
tral energy of the HESS band (lOOGeV-lOTeV). When compar- 
ing the three maps it is immediately apparent that the size of the 
nebula is larger at lower frequency. 

The second feature is that the jet-torus structure, common 
among the observed PWNe in the X-rays and here produced 
by the simulated RMHD evolution in time and space with a 
(Teff > 0.01, is clearly visible also at gamma-ray frequencies. 
One can distinguish the bright features as the central knot and 
the arcs, alt hough the brightness c ontrast is now less than in the 
X-rays (see iDel Zanna et al.ll2006.) . This similarity between the 
gamma-ray and X-ray morphology is due to the fact that the lo- 
cal emitting electrons are the same for both bands. 

An investigation of these bright features in our synthetic 
surface brightness maps at various synchrotron frequencies has 
shown that the central knot and the arcs near the TS appear at 



10 D. Volpi et al.: Non-thermal emission from relativistic MHD simulations of PWNe: from synchrotron to inverse Compton 



IC brightness mop - runA - a,=0.85 - 4GeV 




IC brightness mop - runA - g.^0.85 - 250GeV 




IC brightness mop - runA - 0(,=0.85 - ITeV 




Fig. 5. RunA (2) case. Simulated total surface brightness maps in logarithmic scale normaUzed to the maximum at different IC 
frequency cut-offs: on the left at 4GeV with a maximum 
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all energies and completely dominate the gamma-ray emission. 
Larger structures such as the jets and the torus start to fade away 
from the images at ^ 1 keV and » lOkeV, respectively. Moving 
up in frequency from the synchrotron peak the overall PWN size 
keeps on decreasing reaching a minimum at the cut-off S YN fre- 
quency (O.SGeV in our model), where basically only the elec- 
trons at the TS contribute. On the other hand, the original (radio 
band) dimensions are recovered when IC becomes the dominant 
emission process due to the ubiquity of the radio electrons re- 
sponsible for emission at GeV frequencies. When the wind elec- 
trons begin to dominate the emission, synchrotron burn-off ef- 
fects start to play a role and the volume of the nebula decreases 
again. So all structures, which disappear in the hard X-ray syn- 
chrotron emission, are again visible in gamma-ray images up to 
TeV photon energies. This is due to the fact that the IC emis- 
sion structure does not directly depend on the complex magnetic 
field distribution but only on the spatial behavior of the electron 
distribution function, determined by the local value of p and eoo- 

Present day gamma-ray instruments have insufficient spa- 
tial resolution to distinguish single features in the internal struc- 
ture, so that a comparison between the simulations and the data 
can only be based on the size of the nebula. However, indi- 
cations of asymmetries in PWNe are already seen by HESS 
dMaronian et al. 2005b, 2006a). 

MAGIC has measured an average radius of about 4.5 ly at 
250 GeV and of about 3 ly at 500 GeV dAlbert et alj|2008l) . The 
simulated maps at the corresponding frequencies (we do not dis- 
play here the 500 GeV image) have dimensions comparable to 
the observed values along the y-axis, though along the x-axis we 
find a negligible shrinking of the size with growing frequency. 
This is due to the combined effect of our relatively low mag- 
netic field and the presence of fast flow channels in the equa- 
torial region advecting particles outward fast enough for en- 
ergy losses to be negligible. The average IC fluxes, observed 
by MAGIC, are about half of those obtained in our simulations 
(« 4 X 10"'"erg cm"^ s"'). This corresponds to the excess of the 
synthetic IC emission already discussed when treating the spec- 
tra of Sect. [12l 



5. X and gamma-ray time variability 

In a few PWNe, including the Crab Nebula, bright variable emis- 
sion features are observed at small scales near the termination 
shock and in the torus. In the first case they are called wisps, 
in the second filaments of the torus. The variability time-scale 
is of kilo-seconds in the IR band (Gemini North Telescope, 
Melatos et al.l l2005h . from days up to few months and from 
months u p to a year, at optical and X-ray frequencie s (HST and 
Chandra, IWeisskopf eFallbOOOl: [Hester et alj|2002l) . with time- 
scales that become longer with increasing distances from the TS. 
Quasi-periodically these structures are seen to appear, brighten 
very quickly and eventually fade, showing a characteristic out- 
ward motion. There are no corresponding gamma-ray observa- 
tions because of the insufficient spatial resolution at those fre- 
quencies. 

Different interpretations have been proposed in the liter- 
ature for these ti me- variable features. A kin etic interpreta- 
tion was given by ISpitkovskv & AronsI (|2004|) . who showed 
that if protons are present in the pulsar wind, the vari- 
able structures can arise from ion-cyclotron instability at 
the shock front. Alternative explanations in the context of 
ideal MHD are based on nonlinear Kelvin - Helmholtz insta- 
bilities inside the nebula (Begelman' 4999; 'Bogovalov et dj 
2005; Bucciantini & Del Zanna 2006; Bucciantini et al. 2008|). 
Recently, iBucciantinil (12008) presented simulations with wisps 
near the TS moving outward at velocities of ~ 0.5 c and slower 
filamentary structures in the torus, as observed. The period of the 
variability is of order one year, in agreement with Chandra X- 
ray observations. We are aware that Bucciantini and Komissarov 
are carrying out a detailed study of time variability within MHD 
models, in order to understand its origin, its characteristics and 
possible correlations with other nebular properties. Without go- 
ing into details, we investigate here if the same time-dependent 
flow dynamics responsible for X-ray v ariability can lead to 
fluctuations in the garn ma-ray emission (iDe Jager et al.l 119961 ; 
iLing & Wheaton|[2003h that might be of interest for future in- 
struments (i.e. GLAST). 

We consider runA (2) and compute snapshots of the inte- 
grated spectrum, synchrotron surface brightness map at 1 keV 
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Fig. 6. RunA (2) case. Snapshots from a movie of PWN evolution: output times are selected at 950.6 yr, 955.4 yr and 
960.0 yr. Results from one output time per row are shown. Left panels: plot of the synthetic spectral luminosity in logarith- 
mic scale versus photon frequencies (log). Middle panels: simulated brightness maps (log) at 1 keV normalized to the maxi- 
mum 6.2 X 10"'**erg cm"^ 



sr ' s 'Hz '. Right panels: simulated brightness maps (log) at 250 GeV normalized to the maximum 
1 ur^-i fjjg _^ a[j(j y axes of all surface brightness maps report the distance (in ly) from the central pulsar. 



1.2 X 10" erg cm""^ sr"' s"' Hz 

Notice the inclination with respect to both the plane of the sky and the north direction (respectively 30° and 48° for the Crab Nebula) 

and that only the internal region (within a radius of 3 ly) is displayed. 



(Chandra band), IC surface brightness map at 250 GeV (in the 
middle of the gamma-ray band observed by MAGIC and towards 
the high energy limit of, but still within, the GLAST band). 
The sequence of images lasts ten years, from 950 yr to 960 yr. 



with intervals of 0.2 yr. The full movie is downloadable from 
http ://www. arcetri .astro . it/~delia/crab /runA2.gi^ while a selec- 
tion of snapshots is also shown in Fig.|6] 
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Our simulations confirm the presence of variable wisps near 
the TS and filamentary structures (rings) in the inner part of 
the torus, with characteristic time-scales of about 1-2 years (a s 
observed for energies less than 0.75 Me V, iMuch et aLl 1 1 9951) . 
These features move at speeds 0.3 - 0.5 c (see also the flux ve- 
locity maps in iDel Zanna et al.l2006l) and slow down before fad- 
ing out at larger distances, in agreement with observations. This 
result suggests that variability in the inner nebula can be inter- 
preted as associated to fluid motions occurring in the vicinities 
of the TS (either MHD compressive modes or Kelvin-Helmholtz 
instabilities, the latter arising at the shear flow boundaries). A 
part from these time-varying features, the most prominent struc- 
tures, namely the central knot and the main arc, are instead sta- 
tionary, also as observed. 

The synthetic synchrotron surface brightness maps at X-ray 
frequencies and IC emission maps at 250 GeV photon energies 
show a strongly correlated variability, which is due to motion of 
the common parent electrons. However, as expected, since the IC 
emission is more uniformly distributed in the nebula than high 
energy synchrotron emission by the same electrons, the small- 
scale moving features are less evident compared to X-ray maps. 
The variability of the integrated IC spectrum is correspondingly 
reduced. We would like to recall that in our model at photon en- 
ergies up to 1 TeV, the IC emission is found to be due to wind 
electrons up scattering radio-IR photons, therefore a study of the 
time-variations in this band, still makes sense in spite of the sim- 
plified assumptions underlying our computation of the IC emis- 
sion. We would like to recall that we neglect all propagation 
effects in space and time, approximating the energy density of 
the target photon field as uniform over the nebular volume and 
slowly varying. Both these assumptions are satisfied in the radio- 
IR band. 

For a more quantitative analysis of the variability properties, 
in Fig.|7]we show the synthetic integrated emission as a function 
of time, for a selection of photon energies from X-rays to TeV 
gamma-rays (IkeV, 40MeV, 1 TeV). The strongest variations 
(about a factor of 2) are seen to occur at synchrotron gamma-ray 
frequencies, where the emission is entirely due to the moving 
features close to the TS. On the other hand, the IC time series 
show very limited (about 1%) variations, as anticipated above. 
The results at sy nchrotron frequencies are in agreement with the 
observations by iDe Jager et alJ (Il996h who measured flux vari- 
ations in the 1 - 150MeV band, finding an amplitude of 30% 
in the COMPTEL region (1 - 30MeV) and a factor of 2 in the 
EGRET band (70 - 150 MeV). 

As far as time-scales are concerned, the gamma-ray variabil- 
ity, with a lower limit of orde r a few years in our s imulations, 
also agrees with the results by'Pe Ja ger et alJ (11996) . However, 
the clear relationship between the moving features seen in the 
simulated surface brightness maps of Fig. |6] cannot be easily 
recovered in the luminosity time series, due to spatial integra- 
tion effects. In any case, a combination of XMM/Chandra and 
GLAST observations might confirm the similarities between X 
and gamma-ray variability found in our simulations and conse- 
quently prove their common MHD origin. 

6. Conclusions 



In the present paper the work by lDel Zanna et all (l2006l) on syn- 
chrotron emission from 2-D RMHD simulations of the PWNe 
(and of the Crab Nebula in particular) is extended to the gamma- 
ray frequencies, including the calculation of inverse Compton 
emission self-consistently. For the first time, maps of IC sur- 
face brightness are computed, to compare with existing (e.g. 
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COMPTEL, HEGRA, HESS) and future observations (GLAST) 
of the Crab Nebula. Following recent results in the soft X-rays 
(Bucciantini 2008), a corresponding investigation of time vari- 
ability in the gamma-ray band is also presented. 

In order to com pute the IC emission we include, following 
l Atovan & Aharoni an ( 1996), a homogeneously distributed pop- 
ulation of radio emitting electrons assumed born at the outburst 
of the SN in addition to the wind electrons, continuously accel- 
erated at the TS. In this scenario, the spectral break between IR 
and optical frequencies (^ 10'^ Hz) is intrinsic and due to the su- 
perposition of the two populations, while synchrotron burn-off is 
at the origin of the spectral break in the UV (» 2 x 10'^ Hz). 

From the comparison between our 2-D RMHD simulations 
and the Crab Nebula integrated spectra in Sect. 13.21 we see that 
the value a^ = 0.7, suggested by the optical data of the inner 
pa rts of the nebula (Veron-cet ty & Woltier 1993) and adopted 
bv lAtovan & AharonianI (Il996h . produces an overestimate of the 
synchrotron emission at frequencies above 10'^ Hz. Only with 
the spectral index a,, - 0.85 it is possible to reconcile the 
simulated emission with high energy synchrotron data of the 
Crab Nebula. However the computed IC emission exceeds the 
gamma-ray data by a factor of 2, showing that the steeper parti- 
cle spectrum is not sufficient to remove the discrepancy between 
the results of the simulations and the data. The problem is related 
to the magnetic field in the simulations which is compressed to- 
wards the termination shock and lower on average than the es- 
timates that are found in the literature, and likely than the value 
appropriate for the Crab Nebula. As a consequence, in order to 
reproduce the data we are forced to consider a larger number of 
particles. This seems to solve the problems as far as only the 
synchrotron emission is considered, but its inadequacy is then 
immediately revealed if calculation of IC is added. This sug- 
gests to adopt for the simulations a higher value of the pulsar 
wind magnetization with a word of caution that should be spent 
on the fact that is not clear to us how this would affect the neb- 
ular morphology, which is at present well reproduced. Another 
new result is that the complex nebular dynamics, highlighted by 
our 2-D simulations, gives rise to multi-slopes in the simulated 
synchrotron spectra, which are not expected in 1-D models. 

Further insight can be gained from the gamma-ray mor- 
phology. Our synthetic maps show nebular dimensions and the 
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shrinking of the size with in creasing frequency in the range of 
those observed by MAGIC dAlbert et alj l2008l) only along the 
polar axis, while dimensions do not change appreciably along 
the equator. This is due to the fast flow channels developing in a 
region of low magnetization: these advect the electrons far from 
the termination shock with negligible synchrotron losses. 

The simulated IC surface brightness maps are copies of the 
synchrotron X-ray ones, since they are produced by the same 
parent electrons. However, structures such as the jets and the 
external torus, that disappear in the hard X-rays, survive in the 
entire range of gamma-ray frequencies due to the fact that now 
the emission is not directly affected by the local magnetic field. 

The X and gamma-ray variability of the inner nebula seems 
to have JVIHD origins. We find in our simulations characteristic 
time-scales of about « 1 - 2 yr for both of them and time se- 
ries of the integrated IC luminosity with smaller oscillations and 
slightly longer periods with respect to the synchrotron ones. 

Future developments of the present work will consist firstly 
in treating more accurately the particle evolution in the neb- 
ula, for example by evolving directly the distribution function 
itself along the streamlines. The possibility of a dependence on 
the polar angle along the TS of the injected particle distribu- 
tion function can be thus investigated. On the other hand, in or- 
der to overcome the problem of the magnetic field structure and 
strength, a wider investigation of the parameter space is neces- 
sary (we recall that our results pertain to a single simulation run). 
If this morphology will be confirmed, a problem intrinsic to the 
adopted axisymmetric geometry would be then highlighted and 
3-D simulations should be attempted. 

This analysis could also be extended to other PWNe than the 
Crab Nebula, even in different evolutionary stages (e.g. Vela). 
We believe that the powerful diagnos tic techniques describe d 
here, which complete those presented in lDel Zanna et al.l(l2006l) . 
are sufficiently versatile to be applied to other classes of non- 
thermal emitting sources. 
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